#ifndef ODE_H
#define ODE_H


#include <vector>
#include "linear.h"

using namespace std;




/*vector<MyVec3> ODESolve(vector<MyVec3> state, vector<MyVec3> derivState, float h){
        vector<MyVec3> returnState;
        returnState.push_back(derivState[0] * h);
        returnState.push_back(derivState[1] * h);
        returnState[0] = returnState[0] + state[0];
        returnState[1] = returnState[1] + state[1];
    return returnState;
}*/

MyMat3 ODESolve( MyVec3 L, MyMat3 I, MyMat3 R, float h)
{  
    MyMat3 inverse= MyMat3(6,0,0,0,6,0,0,0,6);
    //inverse.printOut();
    MyVec3 omega= inverse *(L);
    //cout<<L.n[0]<<L.n[1]<<L.n[2]<<endl;
    //cout<<omega.n[0]<<omega.n[1]<<omega.n[2]<<endl;
    MyMat3 omegaSkew=omega.makeSkew();
    //omegaSkew.printOut();
    MyMat3 final= omegaSkew*(R);
    final = final*h;
    //final.printOut();
    return final;
   
}


#endif

